Collective Motion and Phase Transitions of Symmetric Camphor Boats 
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The motion of several self-propelled boats in a narrow channel displays spontaneous pattern 
formation and kinetic phase transitions. In contrast with previous studies on self-propelled particles, 
this model does not require stochastic fluctuations and it is experimentally accessible. By varying 
the viscosity in the system, it is possible to form either a stationary state, correlated or uncorrelated 
oscillations, or unidirectional flow. Here, we describe and analyze these self organized patterns and 
their transitions. 



PACS numbers: 64.60.Cn, 05.65.+b, 45.50.-j 

I. INTRODUCTION 

Groups of motile, interacting bodies are observed ev- 
erywhere in nature, from groups of animals [U-Q to me- 
chanical systems 0, Q . To better understand and unify 
this broad category of collective motion, it is beneficial 
to construct and analyze simpler mathematical models. 
One particularly useful, yet simple model is that of Self- 
Propelled Particles(SPPs) as first described by Vicsek 
and Zafiris [f|. Although the original SPP models are 
very simple, they have displayed a variety of complex be- 
haviors including kinetic phase transitions [6|-|9j and large 
scale pattern formation jlOj . More specific versions have 
also been made, which exhibit similar behaviors 
The benefit of a more specific model is greater applicabil- 
ity to physical systems. Unfortunately, it is often difficult 
to test the results of SPP models because the relevant pa- 
rameters cannot be controlled experimentally. The model 
presented here shares many similarities with SPPs, but 
has the great advantage of being experimentally accessi- 
ble. It displays many of the interesting behaviors seen in 
simpler models as well as new kinds of self organization. 
There is no need to introduce artificial fluctuations, and 
the important parameters are relatively easy to control 
in experiment. 

This model describes the motion of an ensemble of 
symmetric camphor boats(CBs) floating on water. CBs 
are propelled by surface tension gradients generated by 
attached pieces of camphor. They interact, aside from 
collisions, only through their influence on the background 
camphor field. Two-dimensional experiments with sym- 
metric camphor particles have displayed spontaneous 
pattern formation which appears very similar to results 
of 2-D simulations using this model. However, their anal- 
ysis was based on hydrodynamic effects rather than sur- 
face tension. Ensembles of asymmetric CBs have been 
studied experimentally [ill [ljj , but are different in that 
they have a fixed orientation and are driven in a speci- 
fied direction. This paper focuses on the results of 1-D 



simulations and analysis, and some 2-D results are in- 
cluded. In the future we hope to test this model through 
experiment in both one and two dimensions. 

There are three major differences between CBs and ba- 
sic SPPs. First, SPPs always move at a constant speed. 
In this model, and some others [TiTflij ]. boats are subject 
to various forces. However, there is a characteristic speed 
corresponding to a free flowing boat. Second, The direc- 
tion of travel of SPPs is instantaneously set depending 
on neighboring particles. CBs are influenced by collisions 
with other boats and by a surface tension gradient. The 
latter is determined by nearby boats and the history of 
that region. Third, there are no random fluctuations in- 
tentionally added to this system. Most other SPP models 
depend on added noise, and some use it as a key param- 
eter. 

In numerical simulations of the CB model, we observed 
two distinct kinetic phase transitions by varying the vis- 
cosity of the water. One appears as the abrupt formation 
of a stationary, ordered pattern, and the other as a dis- 
continuity in the mean velocity, \{v)\, of the ensemble. 
Aside from these transitions, a variety of collective be- 
haviors were observed such as synchronized formations 
and erratic oscillations. Here we will describe some of 
these behaviors and present quantitative evidence for the 
phase transitions and pattern formation. 



II. MODEL 

The model being considered describes the motion of 
thin plastic disks which float on the surface of water. 
Smaller camphor pellets are attached to the center of 
the underside of the disks. A schematic illustration is 
shown in Fig.([T]). As the camphor dissolves and diffuses 
in the water, it changes the surface tension. The resulting 
tension gradients propel the boats. Such systems have 
often been modeled by a simple set of equations based on 
surface tension and viscosity [13" 20]. This version was 
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a) actual boat 
plastic disk 

camphor disk 



b) 1-D representation 
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c) narrow circular channel 
containing boats 




L=6mm 



FIG. 1. A schematic showing a)an actual boat in 2-D and 
6)the 1-D representation used here, c) shows the periodic 
route in 2-D. 



also used by us in [? ] and these CBs are similar to those 
described in [2l|, but are constructed symmetrically, so 
there is no preferential direction of travel. They move in 
a circular channel with circumference R, which is narrow 
enough to restrict the motion to one dimension. The 1-D 
equations of motion are given by (H}. These equations 
can be readily extended to 2-D by vectorizing them and 
integrating the surface tension around the edge of the 
boat. 

S = --^ + £ h(c(x + L/2)) - 7 ( C (x - L/2))] 
ot z m at m 

(1) 

where m is the mass of the boat, (i is the viscosity con- 
stant of the water and L is the length of the boat. The 
position, x, representing the center of the boat, is defined 
on a periodic domain with period R. The second term 
on the right represents the difference in surface tension 
between the front and back of the boat as a function 
of the camphor concentration given by c(x + L/2) and 
c(x — L/2) respectively. 7(c) is approximated by the sig- 
moidal function in eq|3] 
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Iwater and ^camphor are the surface tension of pure water 
and camphor saturated solution respectively. 

The concentration of camphor molecules on the sur- 
face of the water is constantly changing due to several 
processes. For a system with N boats, it can be approx- 
imated by the following reaction-diffusion equation (J3j) - 

N 



o o2 

" D^-kc + aJ2 F ( x - x i) ( 3 ) 
1=1 



dt dx 2 
F(x) = 1 : for |x| < r , : otherwise 



Here, D is the diffusion constant, k is a constant com- 
bining the effects of evaporation and dissolution, and 
aF(x — Xi) represents the addition of camphor by each 
boat's pellet, which is centered at the point Xi and has 
half-length rrj. To non-dimensionalize the problem, we 
define the following dimensionless quantities. 
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Then the dimensionless parameters of the system become 

fj/ = EIl k' = — r = L3<ylw ~ 7c - > (5) 
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Dropping the ' marks, the non-dimensional equations are 
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F(x) — 1 : for |x| < rp, : otherwise 



To approximate a realistic system, we used parameters 
corresponding to the dimensional values: R — 45.5cm, 
L = 0.6cm, ro = 0.15cm, m = 0.009g, j w = 72g/s 2 , 
7c = 50g/s 2 , D — lem 2 / s (other parameters being var- 
ied). Also, collisions are considered inelastic to match 
the qualitative behavior seen in experiment. 



III. NUMERICAL RESULTS 

In 1-D numerical simulations, three distinct cate- 
gories of behavior were observed. Listed in order of de- 
creasing viscosity, these are: (I)stationary equilibrium, 
(II) oscillation, (III) unidirectional flow. Examples of 
these patterns are shown in the space-time diagrams of 
Fig. @ • The parameter determining the mode of behav- 
ior is the viscosity of the water, /j, in the equations above. 
It also depends on the total density of boats as described 
below. 2-D simulations were also performed, and some 
observations are described, but the behavior of interest 
was more clearly demonstrated in 1-D. 

An order parameter for the transition between station- 
ary and oscillating phases is the e nsemble average of the 
root mean square velocity, yj (v 2 ), where the average in- 
side the root is taken over time. For the transition be- 
tween oscillating and flowing phases, the order parameter 
is the average velocity, \{v)\, which is also averaged over 
both time and the ensemble. The time averages are taken 
over a long interval(Ai > 800), corresponding to hun- 
dreds of oscillation periods and beginning after an initial 
relaxation period. The characteristic oscillation period 
depends on boat density and viscosity. The shortest peri- 
ods were approximately 1 unit of time, but some irregular 
oscillations close to the flow transition took more than 10 
units and required longer simulation intervals. The ini- 
tial positions were semi-randomized to break symmetry. 

We will now qualitatively describe the three different 
types of behavior, enumerated as above, and the tran- 
sitions between them occurring at the critical viscosities 
fx c i and /i c2 . 

(I). In the very high viscosity phase, the system ap- 
proaches a stable, stationary equilibrium in which the 
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FIG. 2. Space-time diagrams showing I) stationary phase, 
II. a) synchronized oscillation, II. b) erratic oscillation, III) uni- 
directional flow (boat density = 0.2). Note the left moving 
dense region in III) representing a jam. hi and h2 are higher 
density versions of II. a and II. b respectively(density = 0.6). 
The lower right figure shows a 2-D stationary state. 



boats are uniformly spaced. 2-D simulations also show 
a stationary, crystal-like pattern which very closely re- 
sembles the experimental observations by Soh, Bishop 
and Grzybowski [l5| . Furthermore, a transition between 
moving and stationary states occurs depending on the 
total density of boats, agreeing with experiment. This 
consistency suggests that the behavior seen in experi- 
ment may be influenced by surface tension rather than 
purely hydrodynamic effects. 

(II). As the viscosity is decreased beyond /z c i, the sta- 
tionary state becomes unstable and the boats begin to 
move. This onset of motion can be quantified by an 
abrupt increase from zero in y/ (v 2 ) as shown in Figj3]a. 
For density greater than about 0.15, the boats oscillate 
and \(v}\ remains almost zero. For lower densities, \(v)\ 
depends highly on the initial conditions, but the change 
in yj (v 2 ) is similar. The 2-D motion depends on the 
configuration of boats and clear oscillation was not seen. 
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FIG. 3. a) yf (v 2 ) in the vicinity of fx c i- b) \{v)\ in the vicinity 
of |Ltc2- Other parameters are: k = 0.072, T = 528, a = 7.2, 
r = 0.25, R = 75, N = 20 



There arc two different oscillation patterns depending on 
viscosity, labeled here as II. o) and II. b). 

II. a) For higher viscosity, the boats oscillate in a very 
synchronized formation as shown in FigJ3] To quantify 
this, we can measure the degree of synchronization be- 
tween boats using the cross correlation of their veloc- 
ities. FigHJa,c shows the correlation vs. the distance 
in numbers of boats. If the density is below about 0.4, 
neighboring boats are anti-synchronized. Higher densi- 
ties show the strongest anti-synchronization at a distance 
of several boats. This distance appears to increase with 
increasing total density. 

II. b) For lower viscosity, the oscillations become more 
irregular and the synchronized behavior vanishes. As 
seen in FigJ3]b,d, the correlation is positive in the vicin- 
ity of the boats and approaches zero for boats further 
away. There is no pattern of synchronization and anti- 
synchronization as seen for type II. a) behavior. 

(III). As viscosity is decreased further, the boats move 
in larger groups and with longer periods between changes 
in direction. There is a critical viscosity, /z C 2, below which 
the boats no longer change direction and the flow be- 
comes unidirectional. The selection of direction is spon- 
taneous and depends on initial conditions. Quantita- 
tively, there is a discontinuity in the net velocity of the 
entire system, \(v)\ as shown in Figj3lb. To test for hys- 
teresis in this discontinuity, the viscosity was gradually 
increased and decreased using a few different values of 
density. For the parameter space tested, hysteresis was 
not observed. This phase can exhibit several different 
kinds of flow including free flow, jammed flow and puls- 
ing flow. Different types of flow have been studied else- 
where [2l|, and they are not distinguished here. 

Fig[5]is a phase diagram illustrating the regions of vis- 
cosity and density for each type of collective behavior. 
A power law fit was made for the density range 0.17 to 
0.38, but it is not clear that the data follow such a rela- 
tion, especially in the case of [l c i- Although the critical 
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FIG. 4. Cross correlation of velocity vs. distance in boat num- 
bers for a) density =0.26, fi = 0.6 b) density=0.26, /i = 4.0 c) 
density=0.6, fi = 0.12 d) density=0.6, fi = 0.68 corresponding 
to behavior shown in Fig[2]II.a, II. 6, hi, h2 respectively 
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FIG. 5. Phase diagram for the types of behavior. The dashed 
lines are a power law fit for densities between 0.17 and 0.38, 
but such a relation is not implied. (same parameters as figO 

values fi c i and [i C 2 are clearly defined, the change be- 
tween patterns II. a and II. & is more gradual. The values 
shown correspond to the point at which the average cor- 
relation between neighboring boats increases to zero. At 
higher density the change was too gradual to assign a 
meaningful transition point, so it is not clearly defined in 
the diagram. The non-stationary phase for density be- 
low about 0.15 depends highly on initial conditions, so 
the distinction between oscillating and flowing phases is 
not shown. However, the transition between stationary 
and moving phases remains even down to the low density 
limit of one boat. 

The transition at /i cl and the change from pattern II. a 
to II. b appear to be independent of the total system size 
as long as the density is kept constant. The value of 
/i C 2 increases very slightly with increasing system size. 
This was tested by varying the total route length between 




10 20 30 40 50 60 70 80 



M 

FIG. 6. Equilibrium velocity for N = 1. 
R = 38 and R = 607. 



IV. ANALYSIS 

The fi c i transition can be analyzed for one boat by 
assuming that it has reached an equilibrium speed, and 
the camphor field is time independent in the frame of the 
boat, c will then satisfy a traveling pulse equation; 

— i>Vc = Ac — kc + aF(x — x$) (7) 

Nagayama, Nakata, Doi and Hayashima have provided a 
detailed analysis of this case for a single camphor pellet 
in 1-D [1 81 ) - They demonstrated a bifurcation between a 
stable stationary solution and a constant-velocity solu- 
tion for a camphor pellet on water. Also, the bifurcation 
is supercritical for small ro, and subcritical for large rn. 
The treatment for a CB is the same, but depends on c 
at the edge of the boat rather than the pellet, so bifur- 
cation type can also depend on the source rate, a. The 
details of this analysis are given in the appendix. The re- 
sults of analytical calculations and simulation are shown 
in Fig|51 This result can be extended to a uniformly 
spaced ensemble of boats by shrinking the periodic do- 
main to include exactly one boat. Unfortunately, it is 
only valid if all boats are moving identically, which does 
not describe the oscillatory behavior observed. In 2-D 
the trajectory of the boat may be irregular as was shown 
in experiment 15], so this technique is not directly ap- 
plicable. 

To explain the oscillatory behavior, there must be an 
effective repulsive force acting on boats that are close 
together. This can be understood by constructing a sim- 
ple scenario. First assume that two boats in the system 
begin moving toward each other, and that their initial 
separation is sufficient that they can approach equilib- 
rium velocity before colliding. Because their speeds are 
nearly equal, and the collisions are modeled as inelastic, 
the boats will nearly come to a stop after colliding. Also, 
assume that the relaxation time of the local camphor 
field is sufficiently short. The camphor field will have 
the form of two neighboring peaks with exponentially de- 
creasing sides. Due to the contribution of each camphor 
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source, the concentration directly between the peaks, cor- 
responding to the closer edges of the boats, will be higher 
than that at the further edges of the boats. The higher 
surface tension on the further edges will provide an effec- 
tive repulsion. The repulsion does not need to have a long 
range, because once a boat is moving in one direction, 
it will tend to continue in that direction. This process 
describes the regular oscillation seen in II. a. However, 
if the velocities of the two boats are significantly differ- 
ent, the final velocity after colliding may be large enough 
that both boats will begin traveling in the same direc- 
tion. The exact speed difference required would depend 
sensitively on the camphor profile and system parame- 
ters. This sensitivity may explain the seemingly chaotic 
oscillation patterns seen in II. b. 

The transition between flow and oscillation cannot be 
treated by such simplifications because the collective be- 
havior is extremely complicated and collision effects are 
significant. Instead, we will outline a reasonable scenario 
based on observations. Consider a single boat in the en- 
semble. The camphor field influencing it can be approxi- 
mated as a combination of the boat's local field generated 
by its camphor pellet and a large-scale field generated by 
the ensemble. The local part typically has the shape of a 
sharp, localized peak, while the large-scale part is much 
smoother and varies over the entire route. Using these 
two parts, the driving force can then be separated into 
local and large-scale parts. The local part typically pro- 
pels the boat in the direction of its velocity, while the 
large-scale part may act in either direction depending on 
the structure of c. We will define /; as the component 
of the force generated by a boat's local field. f g is the 
component generated by the large-scale field. 

Assume that /; is constant, and that f g (x,i) depends 
on the large-scale profile of c(x, t) which may change in 
time. Then the condition /; + f g (x,t) < fxv will cause 



the boat to slow down. Recall that the driving force de- 
pends on the difference in c between the two sides of the 
boat and acts in the direction of decreasing concentra- 
tion. This means that if a boat is moving toward an 
increasing gradient in c, f g will be negative and the boat 
may slow down. Such a condition is often created by a 
boat density peak such as a jam, which generates a peak 
in c. In the flow phase, the boat's momentum is enough 
to overcome the opposing gradient as it approaches the 
peak. If the viscosity is increased, the boat will have a 
lower speed, and may stop and reverse direction before 
reaching the peak. This reversal happens repeatedly for 
viscosity above fi C 2, causing oscillation. Note that the 
boats may pass the peak several times or collide with 
other boats before finally reversing direction. 

In summary, we have theoretically investigated a new 
kind of self-propelled particle and found several distinct 
patterns of collective motion and two kinetic phase tran- 
sitions. The transitions can be quantified by abrupt 
changes in the flow, \(v)\, and the root mean square ve- 
locity, yj (v 2 ) of all the boats in the system. We have 
described the different patterns of self organized behav- 
ior seen in numerical simulations, and outlined a brief 
analysis of the transitions. In future work, we hope to 
test these results through experiment. 
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Appendix A 



x = ±4 



±2r are then 



c(2r ) = B 3 exp (?7_ro) 
c(-2r ) = Ai exp (~i] + r a ) 



(A4) 



The camphor field created by one boat moving at con- 
stant velocity is described by the traveling pulse equa- 
tion. 



Substituting these and the constant velocity, v, into the 
equation of motion gives 



i/Vc = Ac — kc + aF(x — xq) 



(Al) 



The following solution is adapted from Nagayama, 
Nakata, Doi and Hayashima [l8[ It assumes that the 
route length is sufficiently long that the profile of c is 
well localized. This allows us to adopt the boundary 
condition Iimui^oo c(x) — and continuity condition, 
c(x) € C 1 (5ft). Also, we will use the slightly different di- 
mensionless parameters c = — and (3 = (3a to simplify 
the calculation. The solution is 





' Ai exp 


[h+ x \ 


c(x) = < 


A 2 exp 






k B 3 exp 





o = -{iv+r 
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i 



/?B|exp(27 ? _r ) + 1 



(3 A\ exp (-2?7 + ro) + 1 
(A5) 

Note that if v = 0, A\ = B 3 and r; + = — ?7_, so the 
equation of motion is satisfied for any value of /i. We can 
also rearrange the equation to give viscosity as a function 
of equilibrium velocity. 
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se) - T -r 
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where 



Ax 
A 2 
B 2 
B 3 

V = 

v± 



I exp i 



= -2^- 

2k-q 

= eX P (" 

= -f^ eX P 

= n ( ex p ( 

= Vv 2 +4k 
= —v ± rj 



->o) 
- V ~ro) 



exp 



exp 



:vo)) 



(A3) 



(3Bl exp (2ry_r ) + 1 (3A\ exp (-2?7 + ro) + 1 

' A(i) 

This was used to produce the analytical solutions shown 
in FiglBl Inverting this relationship results in a bifurca- 
tion as shown in the figure. By changing (3 in this equa- 
tion, the bifurcation changes between subcritical and su- 
percritical. Recalling the different dimensionless param- 
eters used here, this is equivalent to changing a for the 
parameters used in the main part of the text. 



